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v  A  nine  node  shell  element  is  developed  by  a  new  and  more  efficient  mixed 
formulation.  The  new  shell  element  formulation  is  based  on  the 
Hel linger-Reissner  principle  with  independent  strain  and  the  concept  of  dege¬ 
nerate  solid  shell.  The  new  formulation  is  made  more  efficient  in  terms  of  com 
puting  time  than  the  conventional  mixed  formulation  by  dividing  the  assumed 
strain  fields  into  a  lower  part  and  a  higher  order  part.  Numerical  results 
demonstrate  that  the  present  nine  node  element  is  free  of  locking  even  for  very 
thin  plates  and  shells  and  is  also  kinematically  stable.  In  fact  the  stiffness 
matrix  associated  with  the  higher  order  assumed  strain  plays  the  role  of  a  sta¬ 
bilization  matrix. 
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^  INTRODUCTION 
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As  far  as  the  description  of  geometry  and  kinematics  of  deformation  are 
concerned,  the  advent  of  degenerate  solid  shell  element  [1]  represents  an  impor 
tant  milestone  in  the  development  of  the  finite  element  method  for  analysis  of 
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thin  shell  structures.  The  degenerate  solid  shell  concept  adopts  the  basic 
assumptions  in  shell  theory  which  allows  transverse  shear  deformation  and  an 
isoparametric  representation  in  the  finite  element  approximation.  As  such  it 
can  be  applied  to  finite  element  modeling  of  arbitrary  shell  geometries  without 
resorting  to  a  particular  shell  theory. 

However,  unless  special  care  is  taken,  the  performance  of  a  degenerate 
solid  shell  element  deteriorates  rapidly  as  shell  thickness  becomes  thinner. 

This  undesirable  phenomenon,  called  locking,  reflects  the  inability  of  a  shell 
element  to  represent  zero  inplane  and  transverse  shear  strain  states  without 
disrupting  bending  behavior  [  2] . 

In  an  attempt  to  alleviate  the  effect  of  locking,  reduced/selective 
integration  [3-9]  has  been  widely  used  in  conjunction  with  finite  element  models 
based  on  the  displacement  method.  However,  a  reduced  integration  scheme  does 
not  necessarily  guarantee  a  problem-free  element.  For  example,  2x2  reduced 
integration  rule  applied  to  an  eight  node  displacement  model  element  is  not 
capable  of  eliminating  the  locking  effect  completely.  Moreover,  as  in  the  case 
of  four  node  and  nine  noue  elements  an  excessively  lower  order  of  integration 
introduces  unstable  spurious  kinematic  modes  or  zero  strain  energy  modes  [8-10]. 
In  order  to  suppress  kinematic  modes,  a  stabilization  matrix  can  be  added  to  the 
stiffness  matrix  evaluated  by  a  reduced  integration  scheme  [11,12]. 

Alternatively,  Lee  and  Pian  [2]  showed  that  a  mixed  formulation  based  on 
the  Hel 1 inger-Reissner  principle  or  a  modified  He! 1 inger-Reissner  principle  can 
be  used  to  alleviate  locking.  In  this  formulation,  an  independent  strain  is 
assumed  in  terms  of  parameters  within  an  element  in  addition  to  the  usual 
assumed  displacement.  The  assumed  strain  parameters  are  eliminated  at  element 
level.  Of  course,  for  an  element  with  a  given  number  of  nodes,  the  degree  of 
locking  depends  essentially  on  the  choice  of  assumed  strain.  Moreover,  a  mixed 
formulation  provides  a  rational  mathematical  basis  for  the  reduced  and  selective 


integration  scheme  [2,  13,  14, j.  Following  this  approach,  a  nine  node  shell  ele 
ment  with  specific  assumed  inplane  and  transverse  shear  strains  has  been  deve¬ 
loped  and  applied  to  thin  plates  and  shells  [15,  16].  It  has  been  found  that 
this  nine  node  element  is  free  of  locking  and  has  no  compatible  or  commutable 
spurious  kinematic  modes.  However,  in  the  conventional  mixed  formulation 
used  in  references  15  and  16,  it  is  necessary  to  invert  certain  matrices  in 
order  to  generate  an  element  stiffness  matrix.  Therefore  in  comparison  with  the 
assumed  displacement  model  based  on  the  principle  of  virtual  work,  the  conven¬ 
tional  mixed  formulation  requires  substantially  more  time  to  compute  an  element 
stiffness  matrix  of  the  same  size. 

In  order  to  overcome  this  shortcoming  of  the  conventional  mixed  for¬ 
mulation,  we  develop  in  this  paper  a  nine  node  shell  element  based  on  a  new  and 
more  efficient  mixed  formulation.  This  new  formulation  is  also  based  on  the 
Hellinger-Reissner  principle  with  independent  strain  field.  However  by  dividing 
the  assumed  strain  into  a  lower  order  part  and  a  higher  order  part,  it  is 
possible  to  make  the  new  formulation  much  more  efficient  than  the  conventional 
mixed  formulation.  It  should  be  noted  that  in  reference  17  a  preliminary  study 
on  the  effectiveness  of  the  new  formulation  has  been  made  for  plane  stress 
problem.  In  this  reference,  it  was  demonstrated  that,  for  the  nine  node  plane 
stress  element,  the  new  formulation  requires  less  than  half  of  the  computing 
time  needed  for  the  conventional  mixed  formulation  to  generate  an  element  stiff¬ 
ness  matrix.  In  order  to  more  clearly  demonstrate  the  effectiveness  of  the  new 
formulation  in  comparison  with  the  conventional  mixed  formulation,  we  will  pre¬ 
sent  both  formulations.  The  performance  of  the  new  nine  node  shell  element  will 
be  tested  hy  solving  simple  example  problems  involving  very  thin  plates  and 


GEOMETRY  AND  KINEMATICS 


Figure  1  illustrates  a  portion  of  a  curved  shell  before  and  after  defor¬ 
mation.  For  convenience,  a  local  orthogonal  coordinates  system  is  defined  at  a 
point  on  the  shell  midsurface  in  addition  to  the  global  coordinates  X,  Y  and  2. 
Three  axes  x,  y  and  z  of  the  local  coordinate  system  are  parallel  to  local 
orthogonal  unit  vectors  aj,  32  and  £3  respectively.  The  unit  vectors  aj  and  a£ 
are  tangential  to  the  shell  midsurface  while  £3  is  normal  to  the  surface. 

With  these  coordinate  systems,  the  position  vector  X  of  a  generic  material 
point  P  in  the  undeformed  configuration  can  be  expressed  as 

*sV^3!V?7^  (*> 

where  Xo  is  the  global  position  vector  of  point  0  on  the  shell  midsurface,  t  is 
the  shell  thickness  and  the  nondimensional  coordinate  c  runs  from  -1  to  1.  Note 
that  the  line  connecting  point  0  to  point  P  is  normal  to  the  shell  midsurface.  On 
the  other  hand,  under  the  shell  assumption,  the  displacement  vector  U  of  point  P 
with  respect  to  the  global  coordinate  system  can  be  expressed  as 

“  ■  “0  *  -  7  <*3  -  S3>  <2> 

where  W0  is  the  displacement  vector  of  point  0  in  the  global  coordinate  system 
and  a^  is  the  a3  vector  in  the  deformed  configuration.  For  small  rotation  of 

~3* 

~3  ‘  ~3  =  01  ~1  +  e2  ~2  (3) 

where  rotational  angles  ej  and  63  are  defined  around  the  local  y  axis  and  x  axis 
as  shown  in  Fig.  2.  If  u  represents  the  displacement  vector  of  point  0  on  the 
midsurface  with  components  u,  v,  w  in  the  local  coordinate  system,  the  global 
displacement  vector  U  can  be  expressed  as  follows: 
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where 


U  =  T  u  +  ;  b  9 


(4) 

(4a) 


~  =  ^~1  ~2  ~3^ 

(4a) 

is  the  3x3  transformation  matrix  relating  UQ  to  tj. 

i?  =  7U1 

(4b) 

and 

qt  =  L© j  0^ 

(4c) 

For  a  nine  node  degenerate  solid  shell  element  as  shown  in  Fig.  2,  the 
position  vector  X  of  a  point  in  the  element  can  be  interpolated  from  the  nodal 
values  using  nine-node  Lagrange  shape  functions  Nj(£,n)  such  that 


*  =  j,Ni4+l  XNiM3 


i=l 


i=l 


(5) 


where  X^,  ti  and  a^  are  values  of  XQ,  t  and  a-j  at  node  i  respectively,  and  £  and 
n  are  parent  coordinates.  Similarly,  the  global  displacement  vector  U  is 
assumed  as 


-  *  Wi  +«iwi 


i=l 


i=l 


(6) 


where  u^  and  <3-  are  the  nodal  displacement  and  rotation  vectors  respectively. 

With  the  description  of  geometry  and  displacement  given  above  it  is  easy 
to  establish  a  relationship  between  the  strain  vector  in  ^e  global  coor¬ 
dinate  system  and  the  element  nodal  degrees  of  freedom  vector  ig^.  Neglecting 

p 

higher  order  terms  in  c,  p3  can  be  written  symbolically  as 

f1  =  2q  3e  +  ae  (7) 


r  r  n 

where  B03,  B;.3  are  matrices  relating  p3  to  ge.  Subsequently,  the  strain  vec¬ 
tor  I  defined  with  respect  to  the  local  coordinate  system  is  obtained  by  trans¬ 
forming  the  global  strain  vector  in  Eq.  (7)  as  follows: 


where  Tf  is  the  strain  transform  matrix.  Then  T  can  be  written  symbolically  as 
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f  =  Bg  =  inplane  strain  vector 


(9a) 


s  c  < 


xx 
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r  =  5~<3e  =  bendlng  strain  vector 
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(9b) 
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B  q  =  transverse  shear  strain  vector  (9c) 

~  T  c 


In  the  above  expression  B  B  and  JB  are  matrices  dependent  only  on  e  and  n 
coordinates,  A  similar  expression  holds  for  the  virtual  strain  vectors  sF,  ;6F 
and  6y. 


MIXED  FORMULATION  FOR  SHELL  ELEMENTS 


The  functional  for  the  Hellinger-Reissner  principle  is  expressed  as 


fol lows: 
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E1  C  E)dV 


where,  for  shells, 


1-^xx  ^yy  ^xy  ^yz  ^zx-1 
Independent  local  strain  vector 

^xx  ^yy  ^xy  ^yz  ^zx-^ 

local  strain  vector  derived  from 
displacement  field 
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(10) 


WQ  =  applied  load  term 
Ve  =  volume  of  element 


C  = 


5x5  elastic  coefficient  matrix 


C1  =  3x3  matrix 
C 2  =  2x2  matrix 


The  summation  sign  indicates  summation  or  assembly  over  all  elements. 

In  accordance  with  Eq.  (9),  the  independent  local  strain  vector  E  can  be 
separated  into  inplane,  bending  and  transverse  shear  strains  such  that 


where 


and 
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independent  inplane  strain  vector 


(1 
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XX 


independent  bending  strain  vector 


indenpendent  transverse  shear  strain 
vector. 


(1 


(1 


Moreover,  in  anticipation  of  applying  a  finite  element  approximation,  we  may 
write  dVg  as 

dVe  =  |  J  |  d5  dn  d?  (I 


where  |  J  |  is  the  determinant  of  the  three  dimensional  Jacobian  matrix  J. 


However  |  J|  changes  very  little  along  the  shell  thickness  and  thus  it  is 
possible  to  assume 

dVe  =  )  J  |  ;=0  d£  dn  dc 

Substituting  Eqs.  (9)  and  (11)  into  Eq.  (10)  and  integrating  in  the  c 
direction,  the  functional  becomes 

*R  *  I  [J (£T  £e  £  "  7  £T  £e  £)  dA 

+  /  (i$T  £*  £  "  \  £T  £<  £ )  dA 


+  /  (iTLrVLi)dAl  -  w0 


where 


~e  "  (j  Si  dt 


£<  *  /.  '  Si 


Sy  "  {j  £?  dt 
dA  *1^1  r=0  dE  dn 


The  strain  vectors  e,  ££  and  x  obtained  from  the  assumed  displacement  field  are 
given  in  Eqs.  (9a)  to  (9c). 


(a)  Conventional  Formulation 

In  the  conventional  mixed  formulation,  independent  inplane  strain,  bending 
strain  and  transverse  shear  strain  are  assumed  to  have  polynomial  distributions 
in  terms  of  parent  coordinates  5  and  n.  Symbolically, 

£  “  £e  ~  (14a 
JS-PcU.n)^  H4b 

1  -  P  (5.n)  B  (14c 


where  Pe,  P<>  P^  are  the  shape  function  natrices  of  assumed  strain  fields  and  a 

6^,  6  are  the  assumed  strain  parameters  chosen  independently  for  each  element. 
Substituting  Eq.  (9a)  to  Eq.  (9c)  and  Eq.  (14a)  to  (14c)  into  the  func- 
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where 

Go  =  /  pI  C  B  dA 
~e  1  ~e  ~e  ~e 

(15a 
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no  =--  /  p!  C  Pq  dA 
~e  J  ~e  ~e  ~e 

(15b 
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G  =  /  PT  C  R  dA 

(15c 
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H  =  /  PT  C  P  dA 

~K  **  ~lC 

(15d 
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(15e 
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H  =  /  PT  C  P  dA 
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( 15f 

Setting  6nR  =  0  with  respect  to  a,  BK  and  B  respectively, 

£ 

Go  aQ  -  hq  a  =  o 

~e  ^e  ~e  ~ 

(16a 

•\ 

G  g  -  H  B  =0 

(16b 

G  a  -  H  B  =  0 

*s*y  0  «vy  Aw 

(16c 

These  three  equations  represent  compatibility  in  discretized  form. 

Sol vi ng 

V 

these  equations,  strain  parameter  vectors  a,  8^  and  £  are  expressed 

in  terms  of 

N 

element  nodal  degrees  of  freedom  vector  as  follows: 
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(17a 

► . 

B  =  H"1  G  qQ 

(17b 
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8  »  H"1  G  aa 
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(17c 
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Introducing  a,  8  and  8  into  Eq.  (15),  the  He! 1 inger-Reissner  functional 


can  be  written  as 


"r  -  K^Se  4  se  -  si  5e) 


where 


K  =  GT  H_1  G  +  GT  H-1  G  +  GT  H-1  G 
~e  ~e  ~e  ~e  ~<  ~y 


(18a) 


is  the  element  stiffness  matrix  and 


I  si  3p  -  wf 


(18b) 


with  0  as  the  element  load  vector. 


Since  bending  strain  is  irrelevant  to  locking,  we  may  utilize,  for  the 


finite  element  approximation,  a  modified  Hel  1  inger-Reissner  principle  which  is 
obtained  by  setting  k  =  7  in  Eq.  (13).  A  finite  element  formulation  based  on 


this  modified  Hel 1 inger-Reissner  principle  was  developed  in  reference  15.  For 


the  purpose  of  identification,  the  formulation  based  on  this  modified 


Hel  linger-Reissner  principle  is  called  conventional  formulation  II  while  the 


formulation  described  in  this  section  is  called  conventional  formulation  I. 


(b)  New  Formulation 


The  proposed  new  formulation  is  also  based  on  the  Hellinger-Reissner  prin¬ 
ciple  with  the  functional  given  in  Eq.  (13).  However,  as  suggested  in  reference 
17,  in  the  new  formulation  the  independent  assumed  strains  are  written  as 


E  =  E.  +  Eu 

~  ~L  ~H 


(19a) 


5  =  £L  +  *H 


(19b) 


I  =  XL  +  XH  U9c) 

where  e^,  and  are  independent  strain  vectors  with  lower  order  assumed 
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polynomial  terms  in  c  and  n.  On  the  other  hand,  e^,  and  xH  are  the  higher 
order  strain  vectors  which  contain  higher  order  terms  in  £  and  n.  More 
discussion  on  these  strain  vectors  will  be  given  later.  Inserting  Eqs.  (19a-c) 
into  Eq.  (13),  the  functional  becomes 


where 


rR  =  UI  +  UB  +  US  "  Wo 


Uj  =  I  U(eJ  e  -  \  ej  Cg  £L)dA  +  J  ej  Cg  edA 
-/  4  £e  &LdA  -  7  /  £h  £e  ^HdAl 


(20a) 


K'X' 

Vv' 


few 

-'J 

t 


u8  -  H/ (~i_  £<  £  "  ?  £l  ~ic  £|JdA  +  J^h  ~dA 


■J'-H  ~K  ~LdA  ■  7  ~K  ~HdA^ 


us  =  I  [/ <Xl  £y  1  “  7  II  £y  ^L^dA  +  Jxj  £y  *dA 

■/in  £y  *LdA  ■  7  £y  *HdAJ 


(20b) 


(20c) 


For  the  numerical  integration  of  terms  in  Eqs.  (20a-c),  we  may  use  a  lower  order 
integration  rule  for  the  terms  which  contain  lower  order  assumed  strains  eL,  <L 
and  Xl  respectively.  For  the  remaining  terms  a  higher  order  integraion  rule  is 
used.  Then  the  lower  order  independent  strains  may  be  assumed  in  terms  of 
displacement-dependent  strains  7,  7  and  x  evaluated  at  the  lower  order  integra¬ 


tion  points  as  follows: 


L  ~ 

£L  =  I  N-U.n  )e. 
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c.V.V, 
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(21a) 


wunwrowr 


tt 
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f, 
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\,T 


X 

5b 
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ti 


t: 


£L  -  X  Ni(5,n)£i 


N, 


I  N-U.nW- 

i  =  l  1  1 


(21b) 

(21c) 


where  N,  is  the  shape  function  such  that  *  1  at  the  point  i  of  the  lower 
order  integration  points  and  zero  at  other  points,  F. ,  7.  and  7-  are  evaluated 
at  integration  point  i,  and  N[_  is  the  number  of  the  lower  order  integration 
points.  In  another  word,  we  can  set 


e,  =  e 
~L  ~ 


<  ,  =  < 
~L  ~ 


Xl  sx 


(22a) 

(22b) 

(22c) 


at  the  lower  order  integration  points. 

Introducing  Eqs.  (22a-c)  into  Eqs.  (20a-c)  respectively, 
ui  =  I  It  /  c,  <1A  +  /  el  C  F  dA 


-  ^  ~H  ~L  dA  ‘  7  ~H  ~e  £h  dAl 

ub  =  fL  *[  cK  £L  dA  +  /H  jsJ  F 

‘  ^  ~H  JSl  dA  "  7  ~H  ~k  ~H  dAJ 


(23a) 


dA 


(23b) 


US  =  n?/L  ll  £y  XL  dA  +  /u  ij  Cy  2dA 


*H  £r  iL  dA  •  7  !h  2h  *H  dAJ 


(23c) 


In  Eqs.  (23a-c),  letters  L  and  H  under  the  integral  signs  Indicate  the  lower 
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order  integration  and  the  higher  order  integration  rules  respectively. 

The  higher  order  strains  are  assumed  to  have  higher  order  terms  of  £  and 
n  such  that 

£H=?e{5>Tl)2  <24a 

=5c  (c*n)  £*  <24b 

*H  =  (5‘n)  £  (24c 

where  Pe,  PK  and  PY  are  the  shape  function  matrices  of  the  higher  order  assumed 
strains  and  a,  and  £  are  the  associated  strain  parameters. 

Noting  that  £|_  =  £.  <|_  =  *  and  1|_  =  i  at  lower  order  integration  points, 

and  introducing  Eqs.  (9a-c)  and  (24a-c)  into  Eqs.  (23a)  to  (23c)  we  obtain  the 
following  expressions: 


UI  =  I  (7  Sg  ~!L  ae  +  ST  Ie  3e  -  7  “T  ije  a) 


where 


K 1  -  =  J  bJ  C  B 
~ 1 L  J  ~e  ~e  ~e 


Ga  =  J  Pa  C  B  dA  -  J  P l  C  B  dA 
^6  j ^6  ~e  J . 'vp  ^0  ^0 


He  ■  /„  He  ie  He  dA 


e  ~e  ~e 


where 


ub  =  I  4 3e  £bl  ae  +  £<  ae  -  7  bJ  3*  £k) 

~KBl  -  /  B l  CK  B<  dA 


G  =  J  PT  C  B  dA  -  /  PT  C  B  dA 

~K  4  ^  K  ~K  J  *•'< 


H  *  /  T  C  T  dA 
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and  1Js  =  l  (7  al  J$sL  5e  +  8T  5,  ae  -  7  6T I  6) 


(27) 


where 


Kc.  =  J  BT  C  B  dA 
^-->L  J  ^  ~y  ~T  ~Y 


6  =  /  PT  C  B  dA  -  f  PT  C  B  dA 

~Y  J  ••  ~Y  ~Y  ~Y  '  ,  ~y  ~^y  ~Y 


-Y  ~Y  ~Y 


-Y  "Of  ~Y 


H  =  /  TJ  C  7  dA 

*n/V  /vV 


~Y 


^y 


(27a) 


(27b) 


(27c) 


Substituting  Eqs.  (25),  (26)  and  (27)  into  Eq.  (20),  the  functional  uR  can 
be  expressed  as 

=  I  [7  Se  JSlL  ae  +  aT  Se  -  7  2T  Fe  2 
+  7  3e  ~B|_  ae  +  *1 HK  ae  -  7  gj  B*  g* 

+  7  al  J$sL  ge  +  gT  ^  fle  -  7  ST  ^  g 
"  fll  Op] 


(28) 


Taking  6nR  =  0  with  respect  to  a,  g*  and  £  results  in  three  compatibility  rela¬ 
tions  in  discretized  form  similar  to  Eqs.  (16a-c).  Solving  them. 


a  =  H'1  7  a 
~e  ~b  ^>6 

B  =  H^1  TT,  Qa 
6  =  H^1  U  g 

Substituting  these  equations  into  Eq.  (28)  leads  to 

/I  _T  „  _  _T 


*R  =  l  (7-9e  ~e  3e  ~  3e  2e) 


where  the  element  stiffness  matrix  Ke  is  written  as 

h  a  4  +  *s 


(29a) 

(29b) 

(29c) 
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with 


(31a) 


kl  =  J<IL  +  kBl  ♦  jkSl 


and 


Kc  3  gI  H!1  TT  +  ¥T  H“l  (?  +  (TT  H'1  (T 
~b  ~e  ~e  ~e  ~y  ~y  ~y 


(31b) 


It  should  he  noted  that  the  size  of  ¥e,  ¥y,  lje,  HK  and  ¥y  matrices  in 

the  new  formulation  are  much  smaller  than  Ge,  G*.,  Gy,  He,  H*  and  Hy  in  the  con¬ 
ventional  formulation.  In  addition,  in  order  to  compute "Ge,  H* ,  Tf  an(j  «,  in 

m  /V|  o*  L 

Eqs.  (31a, b),  it  is  necessary  to  evaluate  Be,  B*,  By  matrices  and  |  J  |  at  the 
lower  order  integration  points.  To  save  computing  time  their  values  are  inter¬ 
polated  from  the  values  of  Be,  B* ,  By  and  |  J  |  evaluated  at  the  higher  order 
integration  points  as  follows: 


Se 


Nh  _ 

.i=i  n1  (c,n)  b; 


V 

Bk  *  ,I=i  N.  (c,n)  B; 

v 

By  •  ,I=i  N.  (C,n)  Bj 


(32a] 


(32b! 


(32c! 


V 

I  J  I  »  I  N.  U,n)  |  J  |  1  ( 32d 

1-1  1 

where  Bg,  B^  and  |  J  |  1  are  evaluated  at  the  higher  order  integration  point 
i,  and  is  the  shape  function  such  that  ¥-j  =  1  at  higher  order  integration 
point  i  and  zero  at  the  other  integration  points  and  is  the  number  of  the 
higher  order  integration  points. 

After  assembling  element  stiffness  matrices  and  load  vectors  and  solving 


for  the  nodal  displacement  vector,  the  strain  vectors  are  determined  as 


follows: 


e  =  e,  +7  H  'l  "E  a, 
~  ~l  ~e  ~e  ~e  -'S 


it  s  it,  +  P  H  1  G  q 


Y  s  Yi  +  P  H _1  G  a 
~  ~L  's/Y  '“Y  ~y  ^6 


(33a) 


(33b) 


(33c) 


Then  stress  a  is  determined  from  Eqs,  ( 33a -c )  and  the  stress-strain  relation. 


NINE  NODE  SHELL  ELEMENT 

The  key  to  a  successful  application  of  the  mixed  formulation  is  a  judicous 
choice  of  assumed  strain.  Especially  for  the  purpose  of  eliminating  locking  it 
is  desirable  to  choose  the  simplest  assumed  strain  [2],  However  an  excessively 
simple  form  of  assumed  strain  triggers  spurious  kinematic  modes,  some  of  which 
are  commutable  or  compatible  and  do  not  disappear  even  when  elements  are 
assembled.  This  leads  to  an  unstable  global  model.  Therefore  assumed  strain 
must  be  chosen  such  that  the  element  is  free  of  both  undesirable  locking  and 
spurious  kinematic  modes.  This  often  leads  to  an  element  which  is  not 
invariant.  Invariance  is  not  a  critcal  issue  as  long  as  element  characteristics 
do  not  change  significantly  with  coordinate  systems.  However  in  the  present 
study,  we  enforce  invariance  by  assigning  a  particular  local  coordinate  system 
for  a  given  element  geometry  as  follows: 

For  the  nine  node  shell  element  as  shown  in  Fig.  2,  the  local  coordinate 
system  is  defined  first  by  determining  two  unit  vectors  vj  and  V2  at  £=n=c=0 


point  such  that 


3  X  /  3  X 
-o  /|  ~o 


~1  TC/  H 


(34a) 


3 X  .  3X 
~o  A  ~c 


~2  3n 


/  on 

VK 


(34b) 


The  angle  e0  between  these  two  unit  vectors  is  calculated  from  the  following 


equation: 


(35 


COSB  =  V, 
0  ~1 


Then,  if  60  is  less  than  or  equal  to  90°,  unit  vector  a^  in  the  x  direction  of 
local  coordinate  system  is  chosen  to  be  parallel  with  5  such  that 


3  X. 


a  -  —  / 
~1  '  3£  / 


/ 


H 


(36a 


Otherwise  aj  is  parallel  with  n  such  that 


3Xrt 

a  - 
~1  ~  3n 


3  X 

'-O 


3n 


(36b 


To  determine  the  33  vector  normal  to  the  shell  midsurface,  first  a  unit 


vector  a£  is  defined  as  follows: 


*2 


for  eQ  <_  90°, 
for  eQ  >  90°. 


Then  the  a^  vector  is  obtained  by 

~3  =  ~1  x  ~2  /I  ^1  x  $2  I 


(37a 

(37b 


(38 


and  the  unit  vector  a?  in  the  y  direction  of  local  coordinate  system  is  deter¬ 
mined  as 

~2  =  ~3  x  (39 


Note  that,  while  v^  and  v?  are  determined  at  s=n=5=0  point,  aj,  32  and  £3  can 
be  computed  at  any  point  on  the  shell  midsurface.  Especially  aj,  32  and  53  are 
needed  at  integration  points  to  establish  a  local  orthogonal  coordinate  system. 
Note  that  the  definition  as  described  here  results  in  a  unique  local  coordinate 
system  for  a  given  element  geometry. 


For  the  nine  node  shell  element  based  on  the  conventional  formulation,  we 
assume  the  independent  strain  fields  in  the  local  coordinate  system  as  follows: 


e  =  P  a 
~  ~e  ~ 


where 


P 

~e 


lCnEn  00000000  f 
0000  1  E  n  En  0  0  0  0  0 

0000  OOOOlEnEnO 


and 


2  =  Laj»  a2* 
<  =  P  B 


where 


P  =  P 
~k  ~e 


and 

where 


X  =  L  £ 


P  = 
~y 


8  =  1_8  1 1  8  2> 
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(40) 


(40a) 


(40b) 

(41) 


(41a) 

(41b) 

(42) 


(42a) 


(42b) 


In  Eqs.  (40a)  and  (42a),  f  and  g  represent  higher  order  terms  in  e  and  n  which 
we  will  discuss  shortly. 

It  should  be  pointed  out  that,  for  rectangular  element  shape,  the  conven¬ 
tional  mixed  formulation  element  without  f  and  g  terms  requires  2X2  point  rule 
for  exact  integration  of  element  stiffness  matrix.  Then  the  resulting  element 
is  equivalent  to  the  conventional  assumed  displacement  model  with  2X2  point 
reduced  integration.  This  equivalence  holds  because  the  number  of  integration 
points  is  the  same  as  the  number  of  strain  parameters  in  each  component  of  the 
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assumed  strain  field  and  thus  it  is  possible  to  set  e  =  F ,  <  =  7  and  x  =  I  at 
the  integration  points  [13,14]. 

Note  that  a  nine  node  element  without  f  and  g  terms  in  Eqs.  (40a)  and 
(42a)  will  exhibit  kinematic  modes  or  spurious  zero  strain  energy  modes.  For  a 
flat  square  element  with  sides  along  x  =  ±1  and  y  =  ±1  lines,  these  kinematic 
modes  are  as  follows  [5,15]: 


('  1 

(1) 

u  = 

Cj  x  (1  -  3y2) 

s 

v  = 

-  Cj  y  (1  -  3xZ) 

(43a) 

3 

(2) 

u  = 

c2  (x  +  y  -  3x  y  ) 

(43b) 

(3) 

v  = 

c3  (x  +  y  -  3x  y  ) 

(43c) 

» 

(4) 

w  = 

c4  (x  +  y  -  3x  y  ) 

(43d) 

§ 

(5) 

el = 

c5  x  (1  -  3yz ) 

02  = 

-  c5  y (1  -  3xZ ) 

(43e) 

II 

(6) 

61  ' 

c6  (x  +  y  -  3x  y  ) 

(43f) 

? 

(7) 

e2  . 

c7  (x  +  y  -  3x  y  ) 

(43g) 

> 

where  cj,  ...  , 

C7  are 

arbitrary  constants. 

Among  all  these  kinematic  modes  the 

modes  of  Eqs.  (43a)  and  (43e)  disappear  for  an  assembly  of  only  two  elements. 
However,  all  other  modes  are  compatible  or  commutable  and  they  may  result  in  an 
unstable  finite  element  model  even  after  assembly  of  elements.  Therefore  the 
higher  order  terms  f  and  g  in  Pe,  PK  or  P^  must  be  chosen  appropriately  to  sup¬ 
press  these  undesirable  kinematic  modes.  For  example,  from  Eqs.  (43b)  and  (43c) 


=  C2(2x-6xyZ  ) 


!y  =  C3(2y-6xZy) 


(44) 


Then  the  kinematic  modes  in  Eqs.  (43b)  and  (43c)  will  disappear  if  xyZ  and  xZy 
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terms  are  added  to  the  bilinear  assumed  strains  exx  and  respecti vely. 
Following  the  same  argument,  the  xy  ,  x  y  terms  can  be  added  to  the  bilinear 
assumed  <  and  y  to  suppress  the  spurious  modes  in  Eqs.  (43d),  (43f)  and  (43g). 
Based  on  this  observation,  for  an  element  with  arbitrary  geometry,  we  choose  f 
and  g  as  follows: 

(1)  if  x  is  parallel  with  £  as  in  Eq.  (36a) 


(2)  if  x  is  parallel  with  n  as  in  Eq.  (36b) 
f  =  C2n 

g  =  £n2  (45b) 

Table  1  shows  the  size  of  Ge,  G* ,  &y,  He,  H,. ,  Hy  matrices  which  appear  in 
Eqs.  (15a)  to  (15f)  for  the  conventional  mixed  formulation  I.  These  matrices  are 
evaluated  by  3  X  3  point  integration  rule.  Hue  to  the  need  to  evaluate  these 
matrices  and  invert  He,  and  Hy,  the  conventional  mixed  formulation  I  needs 
more  computing  time  to  generate  an  element  stiffness  matrix  than  the 
corresponding  displacement  model  based  on  the  principle  of  virtual  work.  For 
the  conventional  mixed  formulation  II  based  on  the  modified  Hel 1 inger-Rei ssner 
principle  with  k  =  7,  it  is  not  necessary  to  evaluate  G* ,  and  H"1  .  In  this 
formulation,  the  bending  stiffness  matrix  is  calculated  in  the  same  manner  as  in 
the  case  of  assumed  displacement  model  [15j.  In  the  reference  15,  a  nine  node 
element  has  been  developed  based  on  the  conventional  formulation  II  using  the 
same  assumed  inplane  strain  and  transverse  shear  strain  as  in  Eq.  (40)  and  (42). 
Since  assumed  <  has  nothing  to  do  with  locking,  the  property  of  the  resulting 
element  stiffness  matrix  is  nearly  identical  to  that  of  the  conventional  for¬ 
mulation  I  element.  Therefore  both  conventional  formulation  elements  will  be 
collectively  designated  as  SHEL9. 


On  the  other  hand,  in  the  new  formulation,  we  set  £|_  =  £»  £[_  =  £»  1  =  X  by 
by  adopting  2x2  point  rule  and  a  bilinear  function  for  in  Eqs.  (21a-c)  for 
the  terms  involving  lower  order  integrations.  For  higher  order  integration,  3  X 
3  point  rule  is  used.  Then  the  matrix  in  Eq.  (31)  is  exactly  the  same  as  the 
stiffness  matrix  based  on  the  assumed  displacement  model  with  2x2  point  rule. 
Therefore  it  has  the  same  kinematic  modes  given  in  Eqs.  (43a)  to  (43g).  In 


order  to  eliminate  compatible  or  commutable  modes,  we  choose  the  higher  order 
inplane  strain  as  follows: 


~H  =  ~e  ~ 


where 


V  = 
~e 


f  0 
0  g 
_  0  0  _ 
J  =  Lap  °2  — 1 

For  transverse  shear  strain, 

XH  =  £y  £ 


in  which 


P  = 
~Y 


g  O' 

0  f 


bt  =  Lelt  b2J 


(46) 


(46a) 


(46b) 

(47) 


(47a) 


(47b) 


In  addition,  the  higher  order  is  assumed  to  have  the  same  strain  shape  func¬ 
tion  matrix  as  the  higher  order  assumed  inplane  strain  such  that 


Note  that  the  higher  order  assumed  independent  strain  shape  function  matri¬ 


ces  in  Eqs.  (46a)  and  (47a)  are  taken  from  the  last  two  columns  of  Pe  and  Py 
matrices  in  Eqs.  (40a)  and  (42a)  of  the  conventional  formulation. 

The  size  of  Je,  GK,  Gy,  He,  H*.  Hy  matrices  in  Eqs.  (25)  to  (27)  are  also 
listed  in  table  1.  Due  to  the  small  size  of  these  matrices  as  compared  with 
those  needed  in  the  conventional  formulation,  the  present  new  formulation  will 
require  less  computing  time  for  the  generation  of  an  element  stiffness  matrix 
than  the  conventional  formulation.  In  addition,  to  save  computing  time,  8e,  B* , 

S’  and  |  J  |  at  2X2  points  are  interpolated  from  Be,  B* ,  By  and  j  J  j 
evaluated  at  3X3  points  as  shown  in  Eqs.  (32a-d). 

NUMERICAL  TESTS 

Simple  example  problems  including  both  plates  and  shells  were  solved  to 
test  the  performance  of  the  nine  node  element  based  on  the  new  formulation. 

Plate  bending  problems  are  useful  in  that  the  issue  of  transverse  shear  locking 
can  be  investigated  independently  of  inplane  or  membrane  locking  effect.  For 
the  purpose  of  identification,  the  new  formulation  element  will  be  called 
SHEL9N.  Whenever  possible,  the  performance  of  the  SHEL9N  element  is  compared 
with  the  SHEL9  element  based  on  the  conventional  formulation.  Numerical  results 
are  presented  in  tabular  form  so  that  they  can  be  used  as  a  reference  for  those 

'  's 

who  might  want  to  compare  their  shell  elements  with  the  present  SHEL  9N  element. 
All  examples  were  computed  in  double  precision  arithmatic  on  the  UNIVAC  1100/82 
machine  at  the  University  of  Maryland. 

(a)  Comparison  of  Computing  Time 


turns  out  that  the  computing  time  for  the  new  formulation  element  is  about  50% 
of  the  conventional  formulation  I  element  and  60%  of  the  conventional  for¬ 
mulation  II  element.  Clearly  the  new  formulation  requires  less  computing  time. 


(b )  Simply  Supported  or  Clamped  Square  Plate 

A  quarter  of  a  simply  supported  or  clamped  square  plate  subjected  to  a  uni¬ 
formly  distributed  load  p  is  modeled  by  2  X  2,  3  X  3  and  4X4  meshes.  Both 
evenly  divided  regular  and  distorted  meshes  are  considered  to  investigate  the 
effects  of  element  geometry  and  boundary  conditions  on  the  element.  The  three 
distorted  meshes  are  illustrated  in  Fig.  3.  Elastic  properties  of  the  plate 
are  E  =  107  psi  and  v  =  0.3.  Tables  2  and  3  show  the  nondimensional  maximum 
deflection  at  the  centroid  of  the  plate  for  uniform  regular  mesh.  The  computed 
values  are  normalized  with  respect  to  the  analytical  solutions  at  the  plate 
centroid  obtained  from  thin  plate  theory  [18],  The  results  of  distorted  meshes 
are  given  in  Tables  4  and  5.  The  ratios  of  plate  length  to  thickness  ranges 

2  5 

from  10  to  10  .  The  numerical  values  of  SHEL9  are  taken  from  reference  15. 

The  results  for  3X3  meshes  are  not  reported  in  this  reference.  All  numerical 
results  for  uniform  regular  meshes  are  very  close  to  the  exact  solution  and  are 
insensitive  to  the  wide  range  of  l/t  ratios.  Also  there  is  virtually  no  dif¬ 
ference  between  SHEL9N  and  SHEL9  elements.  For  distorted  meshes,  solutions  of 
both  models  are  slightly  less  accurate  than  those  for  uniform  regular  meshes, 
and  the  performance  of  both  elements  decreases  as  the  plate  becomes  thin,  espe- 
cially  for  the  clamped  plate.  However,  values  of  L/t  over  10  are  beyond  a 
practical  range. 

Table  6  lists  the  nondimensional  maximum  bending  moment  for  the  uniform 
regular  meshes  evaluated  at  the  integration  point  nearest  to  the  centroid  of  the 
plate.  The  computed  values  are  normalized  with  respect  to  the  analytical  solu¬ 
tion  at  the  plate  centroid  obtained  from  thin  plate  theory.  Note  that  for  the 
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conventional  mixed  element,  the  moment  is  computed  at  3  X  3  integration  points 
while,  for  the  present  new  mixed  element,  the  moment  is  computed  at  2  X  2  integra¬ 
tion  points.  Therefore  the  sampling  point  for  the  conventional  mixed  element  is 
closer  to  the  centroid  than  that  for  the  present  new  mixed  model.  Again  the 
reliability  of  both  new  and  conventional  mixed  elements  is  evident. 

(c)  A  Circular  Plate  with  Clamped  Boundary 

A  quarter  of  a  circular  plate  under  uniformly  distributed  load  p  per  unit 
area  is  modeled  by  8-element,  12-element  and  16-element  meshes  as  shown  in  Fig. 


4.  The  radius  R  of  the  circular  plate  is  5  inches  and  elastic  constants  are 

rj 


E  =  10  psi  and  v  =  0.3. 

The  computed  maximum  nondimensional  deflection  w 


64  Dw/pR  at  the  plate 


centroid  is  listed  in  Table  7.  The  analytical  solution  obtained  from  the 


classical  thin  plate  theory  is  wQ  =  1.0.  The  numerical  result  in  Table  7  shows 


that  the  present  SHEL9N  element  is  as  reliable  as  the  SHEL9  element  based  on  the 
conventional  mixed  formulation.  Both  elements  are  free  of  shear  locking  and 
give  accurate  solutions,  especially  for  the  16-element  model.  It  appears  that 
the  SHEL9N  element  solution  is  slightly  more  accurate  than  the  SHEL9  solution. 
However  the  difference  is  very  small. 

(d)  A  Pinched  Cylinder  with  Diaphragmed  Ends 

A  pinched  cylinder  with  diaphragmed  ends  has  been  frequently  used  as  a  deep 
shell  example  because  an  analytical  solution  is  available  for  comparison.  As 
shown  in  Fig.  5,  a  cylindrical  shell  of  length  L  and  radius  R  is  pinched  by  a 
concentrated  load  P  at  two  opposite  points  on  the  circle  of  the  middle  section. 
Due  to  the  symmetry  in  geometry  and  loading,  only  one  eighth  of  the  cylindrical 
shell  is  modeled  by  6  X  4,  7  X  5,  8  X  6,  8  X  6F,  9X7  and  9  X  7F  meshes.  As  an 
example,  the  9X7  and  9  X  7F  meshes  are  depicted  in  Fig.  6.  The  8  X  6F  and  9  X 
7F  meshes  are  obtained  from  7X5  and  8X6  uniform  meshes  respectively  by 
dividing  the  strips  along  the  lines  BC  and  CD  equally.  These  refined  meshes  are 
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used  to  model  the  steep  gradients  of  deflections  and  stresses  in  the  region  of 
the  concentrated  load  point  C.  Material  and  geometric  data  are  E  =  1.05  X  107 
psi,  v  =  0.3,  R  =  4.953  in.,  L  =  2R  and  R/t  =  100,  300,  500. 

The  nondimensi onal  displacements  v7c  =  -wcEt/P  at  the  load  point  C  are 
listed  in  Table  8  for  both  SHEL9N  and  SHEL9  elements  together  with  analytical 
solutions.  Analytical  solutions  were  calculated  by  taking  100  terms  in  each 
direction  for  a  double  Fourier  series  expression  as  given  in  reference  19.  In 
reference  15,  the  9X7  mesh  was  not  tested.  This  table  shows  that  both  mixed 


models  work  very  well  for  this  deep  shell  problem.  Even  though  SHEL9N  has  a 
slight  edge  over  the  SHEL9  model,  the  differences  are  negligible.  For  the  case 
of  R/t  =  100,  the  converged  value  of  is  higher  than  the  analytical  solution 
obtained  by  the  thin  shell  theory  which  neglects  transverse  shear  effect.  In 
Fig.  7,  the  distribution  of  nondimensi onal  deflection  along  the  arc  BC  for  9  X 
7F  mesh  is  compared  for  different  R/t  ratios.  Note  that  curves  for  higher  R/t 
ratios  have  steep  gradients  near  point  C.  This  means  that  more  elements  and 
finer  meshes  are  required  to  model  the  region  near  the  concentrated  load  P  in 


order  to  obtain  accurate  results.  More  numerical  results  are  given  in  reference 

20. 

(e)  A  Pinched  Spherical  Shell 

A  pinched  spherical  shell  shown  in  Fig.  8  serves  as  another  example  of 
deep  shells.  The  spherical  shell  is  a  doubly  curved  shell  in  contrast  to  the 
singly  curved  cylindrical  shell.  The  implementation  of  finite  element  analysis, 
however,  is  relatively  simple  due  to  the  symmetry  in  geometry  and  loading.  In 
addition,  an  analytical  solution  is  available  [21],  Due  to  synsnetry  in  shell 
geometry  and  loading  conditions,  only  a  2°  segment  of  the  upper  half  sphere  is 
modeled  by  four  different  meshes  of  the  SHEL9N  elements.  The  segment  angles  of 
elements  In  different  meshes  are  listed  in  Table  9.  The  region  at  the  loading 
point  is  modeled  by  one  six-node  triangular  mixed  formulation  element  [22]  with 


linear  distribution  of  assumed  inplane  and  transverse  shear  strain  fields. 

Elastic  properties  and  shell  geometry  are  E  =  107  psi,  v  =  0.3,  t=  1  in. 
and  R/t  =  50,  100,  500,  1000.  Table  10  compares  the  nondimensional  normal 
deflection  w  =  Etw/P  at  the  pole.  For  R/t  =  500  and  1000,  the  agreement  bet¬ 
ween  numerical  results  and  analytical  solution  is  very  good.  For  lower  R/t 
ratios,  the  differences  are  slightly  larger.  In  fact,  for  these  ratios  the 
transverse  shear  effect  may  not  be  negligile.  Note  that  Koiter's  analytical 
solution  is  based  on  thin  shell  theory  which  neglects  the  effect  of  transverse 
shear  strain.  Figure  9  represents  the  nondimensional  normal  deflection  of  the 
18-element  mesh.  The  steep  gradient  near  the  pole  is  observed  for  very  high  R/t 
ratios.  More  numerical  results  can  be  found  in  reference  20. 

CONCLUSIONS 


Numerical  test  shows  that,  in  terms  of  computing  time,  the  new  formulation 
is  more  efficient  than  the  conventional  mixed  formulation.  In  fact  the  com¬ 
puting  time  for  the  new  SHEL9N  element  stiffness  matrix  is  about  half  of  that 
for  the  element  based  on  the  conventional  formulation. 

The  SHEL9N  element  is  free  of  locking  and  undesirable  compatible  or  corn- 
mutable  kinematic  modes  as  is  the  SHEL9  element  based  on  a  conventional  mixed 
formulation.  In  fact,  the  stiffness  matrix  associated  with  the  higher  order 
assumed  strain  plays  the  role  of  a  stabilization  matrix. 

The  SHEL9N  element  provides  very  accurate  solutions  for  thin  plates  and 
shells.  Therefore,  we  recommend  it  be  included  in  general  purpose  finite  ele¬ 
ment  programs. 
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Table  4.  Nondimensional  maximum  deflection  of  simply 
supported  square  plate  under  a  uniformly 
distributed  pressure  (distorted  meshes) 
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Table  5.  Nondimensional  maximum  deflection  of  clamped 
square  plate  under  uniformly  distributed 
pressure  (distorted  meshes) 


Nondimensional  bending  moments  Mxx  at  the 
integration  point  nearest  to  the  centroid 
square  plate  under  uniformly  distributed 
pressure  (4X4  regular  mesh) 


L/t 


Type 


Simply 

Supported 


Clamped 


10' 


SHEL9N 

SHEL9 


0.9951 

1.0008 


0.9908 

1.0011 
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SHEL9N 

SHEL9 


0.9951 

1.0008 


0.9907 

1.0012 
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0.9907 

1.0012 
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Table  7.  Nondimensional  maximum  deflection  for  circular 
plate  with  clamped  boundary. 


2R/t 

Type 

Number  of  elements 

8 

12 

16 

102 

SHEL9N 

SHEL9 

1.0023 

1.0024 

1.0014 

1.0018 

1.0013 

1.0015 

103 

SHEL9N 

SHEL9 

0.9982 

0.9931 

0.9948 

0.9943 

0.9979 

0.9974 

104 

SHEL9N 

SHEL9 

0.9805 

0.9632 

0.9822 

0.9809 

0.9918 

0.9906 

SHEL9N 

SHEL9 

0.9790 

0.9546 
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■;*.  Table  9.  Segment  angles  of  elements  for  pinched  spherical  shell 
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Type 

Number  of 
Elements 

Segment  angles  (a$.  ) 
in  degrees  1 

Mesh  1 

9 

2,3,5,7,10,15,20,25 

Mesh  2 

12 

1,1,2,2,4,6,6,8,10,15,15,20 

Mesh  3 

15 

0.5,0.5,1,2,2,4,4,6,8,10,10, 

10,10,12 

Mesh  4 

18 

0.25,0.5,0.5,0.75,1,2,3,4,4, 

6,6,8,8,8,8,10,10,10 
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Fig.  9  Normal  Deflection  along  $  for  the  Spherical  Shell 


